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The Voronoi cell of any atom in a lattice is identical. If atoms are perturbed from their lattice 
coordinates, then the topologies of the Voronoi cells of the atoms will change. We consider the 
distribution of Voronoi cell topologies in two-dimensional perturbed systems. These systems can 
be thought of as simple models of finite-temperature crystals. We give analytical results for the 
distribution of Voronoi topologies of points in two-dimensional Bravais lattices under infinitesimal 
perturbations and present a discussion with numerical results for finite perturbations. 


I. INTRODUCTION 
A. Motivation 

Statistical topology provides a set of tools for studying the statistics of topological properties 
of a system, much in the way that statistical mechanics focuses on geometric quantities such as 
velocities and momenta of large sets of particles. Statistical topology has been applied to many 
classical systems including glasses [T], polymers [2], radial and cellular networks and ideal 

gases [^. In this paper we consider the statistical topology of perturbed lattices. In particular, 
we consider systems in which atoms are initially located at lattice positions and then perturbed 
by random displacements; such perturbations might be associated with thermal vibrations of the 
atoms. In two dimensions, the complete topology of each Voronoi cell is given by its number of 
sides n. What is the distribution of n? We provide exact results for each of the Bravais lattices and 
partial results for the honeycomb structure. Our goal is to provide a more complete understanding 
of realistic systems by providing an analytic point of comparison. 

Statement of Problem. We begin by considering a system of atoms arranged on a two- 
dimensional lattice 


A = {niVi-I-n2V2 I ni,n2 G Z} , (1) 

where Vi and V 2 are linearly independent. A random perturbation of this lattice is obtained by 
displacing each atom from its lattice position by a random variable. Properties of such systems 
are considered in [SIT]. This system can be considered as an approximation for a two-dimensional 
crystalline solid whose atoms are displaced from lattice positions by thermal noise and small strains. 
For a given lattice A, we consider p{n), the probability that a randomly chosen atom will have a 
Voronoi cell with n edges. 

In this paper, we consider perturbations which are chosen from a radially symmetric distribution, 
such that the probability of a particular displacement from a lattice position depends only on its 
magnitude. One particular example of such a model is the classical Einstein solid, in which atoms 
are perturbed from lattice positions using a Gaussian distribution [H]. The Gaussian distribution 
arises from atomic vibrations within a parabolic potential well, where the width of the distribution 
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FIG. 1. Sets of atoms and their Voronoi tessellations: (a) square lattice, (b) perturbed square lattice, and 
(c) random Poisson point process. 


is proportional to the square root of the temperature. We begin by obtaining analytic results for 
infinitesimal perturbations and later consider perturbations of finite amplitude. 


B. Voronoi Tessellations 

We consider Voronoi tessellations in endowed with the standard Euclidean metric d. Let S 
be a discrete, possibly infinite, set of atomic coordinates in Then, 

i?s = {a: € I d{x, s) < d{x, s) for all s' € S} (2) 

is the Voronoi cell associated with an atom s G S. In other words, the Voronoi tessellation is the 
partitioning of the plane into regions such that all points in a region are not closer to any other 
atom than to its own. Note that a point x can be equidistant to multiple atoms and hence belong 
to multiple Voronoi cells. Figure shows Voronoi tessellations for a square lattice, a perturbed 
square lattice, and a set of Poisson-distributed atoms. For our purposes, atoms are said to be in 
general position if no more than three atoms lie on the boundary of an empty circle, i.e., a circle 
that has no other atoms within it. If atoms are in general position, then no point belongs to more 
than three distinct cells, and so no point is adjacent to more than three edges; note that in the 
square lattice this is not the case. In Figs. m and (c), exactly three edges meet at any point. In 
this case, Euler’s formula for planar graphs [S] requires that the average number of edges per cell 
be exactly six. 


C. Classification of Two Dimensional Lattices 

Two-dimensional lattices are traditionally classified into five groups, known as the Bravias lattices 
m and illustrated in Fig.[^ This classification is based on the relative magnitudes and orientations 
of Vi and V2. If Vi and V2 are orthogonal, then the lattice is either square if ||vi|| = ||v2|| or 
rectangular otherwise. If the angle between the vectors is 60° or 120° and ||vi|| = ||v2||, then the 
lattice is hexagonal. For all other angles, the lattice is either rhombic or oblique., depending on 
whether ( 2 v 2 — Vi) • vi = 0. 
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FIG. 2. (a-e) All five two-dimensional Bravais lattices and their Voronoi tessellations, (f) The honey¬ 
comb structure (not a Bravais lattice, since there are atoms not equivalent by translation) and its Voronoi 
tessellation. 


Each lattice has a Voronoi cell which shares the symmetry group of the lattice itself and which 
can be repeated to tile the plane. The square and rectangular lattices have four-sided cells, while 
the hexagonal, rhombic, and oblique lattices each have six-sided cells. The Voronoi tessellation 
of square and rectangular lattices are unstable in the sense that infinitesimal perturbations of the 
atom positions will change the topology of the Voronoi structure. That is because in these lattices, 
vertices in the Voronoi tessellation are equidistant to four atoms. In contrast, the hexagonal, 
rhombic and oblique tessellations are stable, and no more than three edges meet at any vertex. 


D. Delaunay Triangulation 

Every Voronoi tessellation has a dual object called a Delaunay diagram. In a two-dimensional 
Delaunay diagram, edges connect atoms whose Voronoi cells share an edge. If atoms are in general 
position (i.e., no four lie on the boundary of an empty circle), then the Delaunay diagram will be 
a triangulation. Figure [^illustrates a set of atoms, the Voronoi tessellation of the set, and its dual 
Delaunay diagram. 

The Delaunay diagram is a useful tool for computing properties of its dual Voronoi tessellation. 
A Delaunay triangulation maximizes the minimum interior angle of each triangle m- Simple 
mechanisms for finding the Delaunay triangulation, and thus the Voronoi tessellation of a set of 
atoms in general position, are the flip algorithms considered in m- These algorithms begin by 
constructing a triangulation in an arbitrary manner. Next, for each pair of adjacent triangles, the 
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FIG. 3. A set of atoms, the associated Voronoi tessellation (solid lines), and the Delaunay triangulation 
(dashed lines); each Voronoi edge is perpendicular to a Delaunay edge. 

(a) (b) 

FIG. 4. (a) A local non-Delaunay triangulation; (b) the corresponding local Delaunay triangulation. 

edge separating the pair is flipped and it is checked whether this increases the minimal internal 
angle in those triangles. If this is indeed the case, then the edge flip is retained; otherwise it is 
reversed. This ensures that the adjacent triangles are locally Delaunay. This process is repeated 
until no pair of adjacent triangles is non-Delaunay. 

Given the two triangulations for a region bounded by four atoms, picking the one that maximizes 
the minimum internal angle ensures the Delaunay condition is met locally. A simple corollary to 
this theorem is that if a pair of opposite angles adds to greater than 180° on a quadrilateral, then 
the Delaunay triangles split this pair. In Fig. it is clear that (a) has a smaller minimum angle 
than does (b); its is also clear that two opposite angles add to greater than 180° in (a) but not in 
(b). Therefore, Fig. |^b) is Delaunay while |^a) is not. 


II. ANALYTIC SOLUTIONS 
A. Vertex in a Square Lattice 

Before calculating the distribution of Voronoi topologies p(n) for perturbed lattices, we consider 
how the Voronoi tessellation near a topologically unstable vertex (i.e., a vertex with more than 3 
edges) resolves under small perturbations of the atomic coordinates. We focus first on the Voronoi 
tessellation of the square lattice illustrated in Fig.|^a). In the unperturbed case, some points in the 
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FIG. 5. Four atoms and the nearby Voronoi tessellation (shown as solid lines), (a) Before perturbations, 
all internal angles of the Delaunay diagram (Voronoi dual), including a = P, are equal to 7r/2. After 
perturbations, the Voronoi cell and Delaunay diagram change such that either (b) a + /3 < tt and a Voronoi 
edge forms in one direction, or else (c) a + /3 > tt and a Voronoi edge forms in the other direction. 


plane are equidistant to four neighboring atoms and therefore lay on the boundaries of four Voronoi 
cells. After atom positions are randomly perturbed, each unstable Voronoi vertex will resolve into 
two topologically-stable vertices. Such a resolution can happen in two distinct ways, as illustrated 
in Figs. I^b) and (c). Given the symmetry of the problem, it seems intuitive that each resolution 
should occur with equal probability. Here we make this intuition precise. 

As noted above, the Delaunay diagram can be used to study properties of its dual, the Voronoi 
tessellation. In the unperturbed case illustrated in Fig. |^a), all internal angles of the Delaunay 
diagram are tt/2, including the pair a, j3. After a random perturbation, a+/3 will almost surely either 
decrease or increase, determining whether the unstable vertex resolves in the manner of Fig. ib) 
or[^c), respectively. Determining the probability of obtaining these configurations thus requires 
calculating the probability with which a + j5 increases or decreases under a random perturbation. 

In the remainder of the paper we denote a spatial configuration of n atoms in as 


X = 



X2 

. 2/2 _ 


.Vn\J 


(3) 


such that the initial, unperturbed configuration of four atoms on a square lattice is 


X = 


(4) 


We use (/>(X) to denote the sum of the angles a and j3 illustrated in Fig. Ufa). We apply a random 
perturbation to (j), namely ^ G where ^ is chosen from any radially symmetric probability 

distribution. As (f) changes smoothly with atomic coordinates, the small perturbation p can be 
linearized as Vt/) • such that (j){X + ^) « 4’i^) + '^4’ ‘ ?• 

In general, the set of perturbations p that satisfy ^ = 0 can be represented as a hyperplane 
in passing through the origin. For perturbations p on one side of this hyperplane we have 
• ^ > 0, while for those on the other side we have Vt/) • ^ < 0. For the configuration to resolve 
as shown in Fig. ib), it must be that ^(X + ^) > tt or, equivalently, V(/) • ^ > 0. If Vt/f • ^ < 0, 
the initial configuration will resolve as shown in Fig. ic). Since the distribution of perturbations 
is symmetric about the origin, the integral of the probability measure over the region for which 
V(^-^>0is0.5; this is also the integral over the region for which ^ < 0. This demonstrates 
that the topological resolutions illustrated in Figs, [^b) and (c) indeed occur with equal probability. 
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(pi — Oil + /?1 

</'2 = ce 2 + /32 

(a) 


V01 ■ ^ > 0 
V (/)2 ■ C > 0 
(b) 





C > 0 

V</)2 ■ ^ < 0 

(c) 



V</)i ■ ^ < 0 

V </)2 • C < 0 
(d) 


FIG. 6. (a) Prior to perturbation, all internal angles of the Delaunay diagram are tt/2, and the two vertices 
of the Voronoi tessellation are unstable. After perturbation, each vertex will resolve in a manner that 
depends on Oi + /3i. (b-d) show how the topology of the Delaunay diagram resolves depending on the signs 
of Vt/ii • 5 and V(^2 • C- Flipping the signs in (c) leads to the remaining resolution, which is identical to (c) 
up to rotation by 180°. 


In this example, the probabilities depend on the fractional volume of space on either side of a 
single hyperplane. We next consider probabilities of more complex configurations, the calculations of 
which require consideration of fractional volumes associated with multiple intersecting half-spaces. 


B. Nearby Vertices in a Square Lattice 

Another example illustrates the complexities that can arise from correlations between neighboring 
unstable vertices. Figure illustrates two nearby, unstable vertices in a Voronoi tessellation of a 
square lattice. The manner in which these resolve is determined by the motions of seven nearby 
atoms, and how the labeled angles at and /3i change under perturbations. 

Before perturbation ai = /3i = tt/2 for all i, and the two vertices of the Voronoi tessellation 
are unstable. After perturbation, if -I- A > tt a new edge will form between the two atoms, as 
illustrated in Fig. |^c). As shown above, the probability of -b ft > tt is 0.5. What is the joint 
probability that both ai + j3i and 02 + ft will increase under a random perturbation? 

As before, we denote a configuration of 7 atoms as a vector X G where n = 7; we let 
(j)i(X) = ai + Pi- Each of ft = 0 determines a hyperplane through the origin in K^", and 
the two intersecting hyperplanes divide into four regions, each corresponding to one of four 
possible topologies that result from the infinitesimal perturbations of the atoms, as illustrated in 
Fig.^b). If we let the initial, unperturbed configuration space of atoms be: 



1 0 
1 ’ 0 


-1 

-1 



V = 


(5) 
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FIG. 7. Schematic of divided into regions in which the topology resolves in distinct manners. Each 
W4>i divides R^" into two regions, one on which • ^ < 0, and the other on which V(?ii • ^ > 0, where ^ is 
a point in R^". 


then the gradients 'S/4>i are: 


V</.i 


V (^2 


'O' 


'I' 
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'O' 

1 

1 
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-1 
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5 
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_ 1 

5 

- 1 
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and 


0 

-1 


( 6 ) 

(7) 


The two hyperplanes defined by these gradient vectors intersect at an angle 9 = arccos(—1/4) « 
104.48°; this angle does not depend on the choice of coordinates. 

Since the perturbations are chosen from a probability distribution that is radially symmetric 
about the origin, the probabilities of the configuration resolving into each of the four possibilities 
are equal to the fractional volumes of the four regions separated by the intersecting planes (see 
Fig.0. These fractional volumes can be determined through consideration of the angles at which 
these two planes intersect. 

In the case considered here, the angle between Vi/i and V ^2 is 0 = arccos(—1/4), and so the 
angle subtended by the region where Vi/i • ^ > 0 and — 0. The fractional volume of 

this region is then (tt — 9) /27r « 20.98%. The fractional volume of the region where V(/i • ^ < 0 and 
V (/2 • ^ < 0 is the same (see Fig. |^. Hence, the probability of the initial configuration resolving 
into the configuration shown in Fig. [^b) or (d) is (tt — 0 )/ 27 r for each. By the same reasoning, 
the probability of the initial configuration resolving into the configuration shown in Fig. [^c) is 
9/2tt « 29.02%. 


C. Square Lattice 


The above approach allows us to compute p{n) - the probability that the Voronoi cell of a 
randomly chosen atom has n edges - in a perturbed square lattice. Figure illustrates a central 










































atom in a square lattice and its eight nearest neighbors. The resolution of all four unstable vertices 
in the Voronoi tessellation can be determined by considering how the sums of angles (pi = ai + Pi 
change under perturbation. Following the approach above, the resolution of the entire configuration 
under a perturbation ^ is determined by the signs of 


The resolution of each unstable vertex can be described by a pair of half-spaces in with 
n = 9; in one half space ^ > 0, and in the other Vpi ■ ^ < 0. The way in which all four 
unstable vertices resolve is then described by the intersection of four of these half-spaces - one for 
each pair. Under a random perturbation the Voronoi cell of a central atom gains one edge for 
each i for which Vpi ■ ^ > 0. The probability that the Voronoi cell of an atom in a square lattice 
resolves with 8 edges is the fractional size of the region in in which • ^ > 0 for all i. The 
probability that a Voronoi cell resolves with 4 edges is equal to the fractional size of the region in 
in which Vpi • ^ < 0 for all i. 


If the initial, unperturbed configuration of the atoms is 


V = 
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( 8 ) 



FIG. 8. (a) An atom and its eight nearest neighbors in a square lattice; (b) this atom under perturbation, 
where the resolution of the Voronoi tessellation is determined by the sums at -\- Pi formed by neighboring 
atoms, as illustrated here and described in the text. 
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the gradients Vfjji 

V,^i= ( 
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V02 = 
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( 10 ) 


v^3 = 
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V(^4 = 
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( 12 ) 


By consideration of their dot products, it is straightforward to see that V(()i is orthogonal to V ^2 
and V(^ 4 , though not to Wcj)^. In fact, each W(j)i is orthogonal to the W(j)i associated with its two 
nearest neighbors. For any two orthogonal vectors wi, W 2 and a random perturbation ^ chosen from 
a radially symmetric distribution, P[wi > 0, W 2 - ^ > 0] = F’[wi > 0 ]P[w 2 > 0]- This follows 

from the fact that a radially symmetric distribution is invariant under rotation. Therefore, we can 
compute the joint probability ^ > 0, V (/)3 • ^ > 0] independently of the joint probability 

P[V02 • C > 0, V (/)4 • ^ >0], and so the probability of a Voronoi cell resolving with 8 edges is given 
by: 


p{8) = P[V(^, • C > 0] 

= P[V(^i-^>0,V<()3-C>0] 


P[V<(.2-e>0,V<()4-e>0]. 


(13) 


Since we have shown (Section |ll B ) that P[V(j)i ■ ^ > 0, V (^3 • ^ > 0] = P[V (()2 • ^ > 0,V(j)4 ■ ^ > 
0] = (tt — 6)/2tt, the joint probability of V(()i • ^ > 0 for all i is [(tt — 0)/27r]^ « 4.40%. Because the 
probability distribution of the perturbations is radially symmetric, this is also the probability p(4). 
This is consistent with numerical results reported below in Section |III[ 

We determine p{5) and p(7) in a similar manner. For the Voronoi cell to resolve with five edges, 
one vertex must resolve to gain an edge and all remaining vertices must resolve so that no other 
edge is gained. This can occur in four ways, one for each vertex that can gain an edge; the sum of 
the probabilities of these four events is: 


p{5) = • C > 0, V<(.3 • ^ < 0] • P[V<(.2 • e < 0, V((.4 • C < 0] + 

• C < 0, V</>3 • ^ > 0] • P[V</.2 • C < 0, V(/.4 • C < 0] + 

P[V<(.1 • C < 0, V</>3 • ^ < 0] • P[V</>2 • C > 0, V(/.4 • C < 0] + 

P[V<(.1 • C < 0, V</>3 • ^ < 0] • P[V</.2 • C < 0, V(/.4 • e > 0] (14) 

= « 24.35%. (15) 

Noting that the probability of two opposite sides both gaining an edge is the same as both not 
gaining an edge, clearly p{5) = p{7). We can calculate p{6) by noting that we have already 
determined every other case (since it must resolve so that the central atom has 4 to 8 edges); i.e., 

p(6) = 1 - p(4) - p(5) - p(7) - p{8) = i - f + i (f)' « 42.49%. 
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D. Rectangular Lattice 


Our approach for the square lattice naturally leads to results for rectangular lattices. We can 
describe the aspect ratio of a rectangular lattice as y = |lvi|l/||v 2 ||, where Vi and V 2 are basis 
vectors for the lattice as described in Section [TCl Following the conventions and definitions for the 
square lattice, if the initial positions of the atoms are given by the configuration 
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then the gradients X(j>i are 
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(16) 

(17) 

(18) 

(19) 

( 20 ) 


Under a random perturbation the rectangular Voronoi cell of a central atom will gain an edge 
for each i for which ^ > 0, in a manner similar to that illustrated in Fig. Determining 

the distribution of Voronoi topologies p{n) of a perturbed rectangular lattice requires computing 
the fractional size of the region in which ^ have particular signs. In the general case, this 

fractional volume is equivalent to the solid angle subtended by the region, divided by the surface 
area of the unit n-sphere. 

No closed-form expression exists to calculate the solid angle bounded by n linearly-independent 
vectors when n > 3 m- Nevertheless, several approaches are available to numerically compute 
this solid angle m- In this paper, we compute the probabilities p(n) for each n by numerically 
integrating the Gaussian function over the intersection of the four half spaces. Although any radially 
symmetric probability function could have been used in the integration, the Gaussian was selected 
due to its well-understood functional form. Figurej^shows p{n) obtained from numerical integration 
for aspect ratios 0 < y < 10. Note that p{n) is the same for y as for 1/y. 

It is interesting to inquire how p(n) behaves as y approaches 0 or oo. Limit analysis of the 
integral equations, combined with numerical integration, shows that p(4) = 5.46%, p(5) = 20.09%, 
p{6) = 48.90%, p{7) = 20.09%, and p{8) = 5.46% as j/ —>■ oo or y —0 (see Fig. |^. 


E. Vertex in a Honeycomb Structure 

We next focus on the honeycomb structure shown in Fig. [^f); we note that this is not a Bravais 
lattice, since there are atoms not equivalent by translation. In the center of the hexagonal arrange¬ 
ment of atoms is an unstable Voronoi vertex which is adjacent to six edges. With probability one. 
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FIG. 9. Numerical integration results for an infinitesimal perturbation as a function of rectangular lattice 
aspect ratio, y. When the aspect ratio is one, the result is the square lattice. Large circles indicate p(n) 
for aspect ratio 1, and also for asymptotic values as the aspect ratio approaches infinity. 


the original unstable vertex resolves into three stable vertices, each of which is adjacent to three 
edges. Such a resolution can occur in one of three ways, as illustrated in Figs. 10 |a-c). 


To compute the probabilities of each of these resolutions, we consider the initial configuration of 
atoms 


X = 


and the corresponding gradients 
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( 21 ) 
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( 23 ) 

( 24 ) 


where (j)i = ai + j3i, for a resolution of the initial configuration as shown in Fig. if)- Each of 
the resolutions shown in Figs. 10 a-c) has a triplet of gradient vectors associated with it that 
determine the probability of the initial configuration resolving as shown. The Voronoi topology 
resulting from a perturbation can be determined through examination of the Delaunay diagrams. 
Each topological resolution corresponds to a distinct Delaunay triangulation, the topology of which 
depends on satisfying a set of inequalities in terms of internal angles formed by the six neighboring 
atoms. 

For example, in order for the unstable vertex to resolve in the manner illustrated in Fig. &)> 
it is necessary that > tt for all i; this occurs when W(j)i ■ ^ >0 for all i. The probability of 

obtaining this resolution is equal to the fractional region of space of three intersecting half spaces. 
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FIG. 10. (a-c) Resolutions of unstable vertex in the honeycomb structure; (d-f) corresponding Delaunay 
diagrams. 


each described by • ^ > 0. We can compute this fractional region by considering fl, the solid 
angle subtended by this region of space. For the configuration described by the Voronoi tessellation 


in Fig. 10 |a) and the Delaunay the in Fig. 10 d), we compute = 2 tan ^ (-\/6 — 2). We then have: 


c|v«...{>oi = F = 


2 tan ^ (Vo — 2) 




0.067, 


(25) 


since the fractional volume of a three dimensional wedge is its solid angle divided by the surface 
area of a three-dimensional sphere. Accounting for symmetries, Fig. [^a) can be obtained in six 
equivalent ways (obtained by 60° rotations), such that the probability of obtaining this topology is 
6 £. > 0] « 40.34%. Similar calculations show that the probability of obtaining the topology 

(^ 


in Fig. 


10 F) is ^ tan 


2 V6-3 


34.66% and in Fig. 


10 


c) is 1/4. 


F. Hexagonal, Oblique, and Rhombic Lattices 


For infinitesimal perturbations, the hexagonal, oblique and rhombic lattices are all stable, since 
for any empty quadrilateral, opposite angles never sum to tt; the Delaunay condition (Section ID I 
shows that these are all stable triangulations. With probability 1, an infinitesimal perturbation 
cannot change this; in the limit that the magnitude of the random perturbation goes to zero, there 
will be a point at which the linearization |V(/) • ^| < 7 where 7 is the difference from the original 
sum of opposite angles and tt. Therefore, in this limit, these structures are stable. As an example, 
consider the hexagonal lattice from Fig. if)> where any empty quadrilateral has a pair of opposite 
angles that sum to 27r/3 and thus are stable for all small perturbations. 

Table ^ summarizes the analytic results from Section ^ as well as the known Poisson-Voronoi 
distribution. 
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P(3) 

p(4) 

p(5) 

p(6) 

P(7) 

P(8) 

p(9) 

Square (exact) 

0 

[i^-e)/2nr 

6{7r-e) 

TT^ 

1 G , 3 ( eV^ 

0(7r-e) 

TT^ 

[(^ - e)/2nf 

0 

Square (approx.) 

0 

0.0440 

0.2435 

0.4249 

0.2435 

0.0440 

0 

Rectangular (3) 

0 

0.0499 

0.2198 

0.4605 

0.2197 

0.0499 

0 

Rectangular (10) 

0 

0.0539 

0.2033 

0.4856 

0.2033 

0.0539 

0 

Rectangular (oo) 

0 

0.0546 

0.2009 

0.4888 

0.2009 

0.0546 

0 

Hexagonal 

0 

0 

0 

1 

0 

0 

0 

Rhombic 

0 

0 

0 

1 

0 

0 

0 

Oblique 

0 

0 

0 

1 

0 

0 

0 

Poisson 

0.0112 

0.1068 

0.2595 

0.2947 

0.1987 

0.0897 

0.0295 


TABLE I. Summary of results for the distribution of Voronoi topologies under infinitesimal perturbation for 
each Bravais lattice and for Poisson distributed points M- As explained in Section[lIB[ 6 = arccos(—1/4). 


III. SIMULATIONS 


While we reported above on the effects of infinitesimal perturbations on the topological distri¬ 
butions of Voronoi cells, we now examine the impact of the perturbation amplitude on the same 
distributions. For this purpose, we displace each atom i from its initial lattice position Xi by a 
two-dimensional independent random variable Axi chosen from a Gaussian distribution centered 
at the origin, and described by the probability density function 

/(Ax.) = (26) 

the new position of each atom after the perturbation is x' = Xi -I- Ax^. 

We consider how the distribution of Voronoi topologies p{n) in the perturbed system varies 
with the width of the Gaussian a. When a is small, the results obtained in Section should 
be closely reproduced. When a is large, p{n) approaches the Poisson-Voronoi distribution, as the 
distribution of atoms within any region approaches a Poisson distribution. The Poisson-Voronoi 
distribution is known analytically in integral form |14j and has been evaluation numerically and 
through simulations m- 


A. Simulation Design 


All simulations were designed in Matlab using the standard Voronoi package and random number 
generator. A grid of atoms for a lattice structure was initialized and each atom was displaced from its 
lattice position by choosing random displacements in x and y according to a Gaussian distribution 
with mean p = 0 and standard deviation a. The standard deviation a was measured in units of 
the square root of the unit cell area: v/||vi x V 2 ||. 

Each simulation considered a system of at least one million atoms; periodic boundary conditions 
were used to minimize finite-size effects. The simulations recorded the fraction of Voronoi cells with 
a particular topology. 
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FIG. 11. Frequencies of n-sided Voronoi cells for 0 < ct < 1 for the square lattice; large circles indicate 
asymptotic values. 


B. Square Lattice 


Figure 11 shows the distribution of Voronoi topologies p{n) of a perturbed square lattice for 
0 < CT < 1, obtained through numerical simulation as described above. Simulation data for small 
a are consistent with the analytical results reported in Section |II C| The symmetry of the analytic 
result is quickly broken as a increases. Whereas p(5) = p{7) and p(4) = p{8) in the limit as a 
approaches 0, for finite a, this is not the case. Moreover, where p{n) = 0 for n = 3 and n > 8 in the 
small a limit, these numbers become positive for finite a. As a increases, p{n) evolves monotonically 
toward the Poisson-Voronoi distribution except for p(4) which has an inflection point near a — 0.3. 
The fraction of Voronoi cells with n > 8 edges increases rapidly between a = 0.2 and a = 0.4, and 
then levels off by ct = 0.5. By ct = 1, it is difficult to distinguish the distribution p{n) from the 
Poisson-Voronoi distribution. 


C. Rectangular Lattice 


Figure 


12 shows the distribution of Voronoi topologies p{n) of rectangular lattices of aspect ratio 


three and ten respectively for 0 < ct < 1. For low ct, the numerical results are consistent with the 
analytic results obtained for infinitesimal perturbations. Qualitatively, the Voronoi cells with 4, 6, 
and 8 edges are favored more heavily than they are in the square lattice. Voronoi cells with 5 and 
7 edges become more frequent until an abrupt change around ct = 0.2 for aspect ratio three and 
CT = 0.1 for aspect ratio ten. While p{n) is monotonic for most n in the square lattice, p{n) is not 
monotonic for most n in the rectangle lattice at aspect ratio ten. As the aspect ratio increases, the 
rectangle lattices show more significant inflections between the low and high aspect ratio limits. 
For example, there is a sudden peak of three-sided Voronoi cells near ct = 0.5 that disappears at 
higher ct. The convergence to the Poisson-Voronoi distribution is also slower as the aspect ratio 
increases. 
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FIG. 12. Frequencies of n-sided Voronoi cells for 0 < cr < 1 for rectangular lattices of aspect ratio 3 and 10 
respectively. 
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FIG. 13. Frequencies of n-sided Voronoi cells for 0 < n < 1 for the hexagonal lattice; large circles indicate 
asymptotic values. 


D. Hexagonal Lattice 


Figure 13 shows the distribution of Voronoi topologies p(ri) for the hexagonal lattice, obtained 
through numerical simulation, as a function of a. The distribution p(n) is relatively insensitive to 
changes of a at both low and high a limits. The switch between these limits is abrupt; it occurs 
largely within the limits 0.1 < a < 0.3. As a result, the distribution p(n) varies monotonically for 
most n between the small and large a limits. 
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FIG. 14. Frequencies of n-sided Voronoi cells for 0 < ct < 1 for the honeycomb structure; large circles 
indicate asymptotic values. 


E. Honeycomb Structure 

Figure [l^ shows the distribution of Voronoi topologies of the honeycomb structure for 0 < cr < 1 
as obtained through numerical simulation. The variation of p{n) with cr is “gentle” relative to 
the lattices discussed above. This might result from the similarity between the low- and high-cr 
frequencies p(n). In the honeycomb structure, even as a approaches 0 there is a large diversity of 
topological types, with non-zero frequencies for n = 3 and 9 < n < 12; in contrast, for other lattices 
these frequencies are all zero as a approaches 0. 


IV. CONCLUSIONS 

In this paper, we explored how perturbations affect the distribution of Voronoi cell topologies for 
the two-dimensional Bravais lattices and the honeycomb structure. Such results are important for 
describing crystalline materials at finite temperature. Of particular interest are those structures for 
which the zero-temperature Voronoi topology distributions are unstable with respect to infinitesimal 
perturbations. In applications to crystalline materials, it is the small-perturbation limit that is 
of most interest at temperatures well below the melting point. Even approaching the melting 
temperature, a rarely exceeds 0.05 — 0.15 m- 

Voronoi cell descriptions of the local structure have been growing in importance for material 
characterization, both for atomic-scale simulations and in experiment. Although here we focus 
on two-dimensional lattices, the same methods can be applied in three dimensions, where most 
materials and structures live. Indeed a similar approach has been used in recent work [17] to 
characterize disorder in three-dimensional atomic systems. 

While other papers have focused on topological distributions for Poisson distributed points, our 
paper considers the topology in systems which have much more order, but are not strictly periodic. 
The most significant achievement in this paper is the summary of exact results for the distribution of 
Voronoi topologies of all two-dimensional Bravais lattices under random infinitesimal perturbations 
chosen from any radially symmetric distribution (see Table . Because we present numerical data 
for finite perturbations, these results also show how the topology transitions from that of a perfect 
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lattice to a perfectly disordered (Poisson) system. While the present study focuses on uncorrelated 
perturbations in the Einstein-sense, it would be of interest to see how accurate this description is 
for realistic atomic systems in which vibrations are correlated, as described by the phonon spectra. 
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